# Full text of "DIRECT CALCULATION OF SUBLIMATION ABLATION"

## See other formats

> NASA TECHNICAL TRANSLATION NASA TT F-l4,529 DIRECT CALCULATION OF SUBLIMATION ABLATION H. Kubota Translation of: "Shoka abureshon no chokusetsu kalho", Koku Uchu Gijutsu Kenkyusho Hokoku, [Technical Report of National Aerospace Laboratory] , To K I O^ No. 239, June, 1971, PP- 1-25 SlSftfon ""fee? «u.. ^"^ ^° " "fo. G3/33 N72-32936 Unclas U1014 ;♦' fjiSSzi?^ kO> >S J^nZ'S?: NATIONAL AERONAUTICS AND SPACE ADMINISTRATION WASHINGTON, D.C. 205^6 AUGUST 1972 Reproduced by' SPnngf,e|JvA%i?,'«» DIRECT CALCULATION OF SUBLIMATING ABLATION Hlrotoshi Kubota Abstract. A direct and rapid calculatlonal . method of sublimating ablation for a hypersonic '■ heat-shielding device is presented. This is based on a numerical analysis of the Hartree- Womersley method which reduces parabolic partial differential equation to a corresponding ordinary one by means of a kind of implicit. 1. Introduction When a spacecraft reenters the atmosphere, powerful shock waves /L*** are produced around the craft, and the air currents behind the shock|waves assume an extremely high temperature. It has been made clear in recent research that an effective way to protect the craft from this sort of aerodynamic heating is to adopt mass transfer cooling, in which a relatively small quantity of cooling gas or cooling liquid is blown off into the external fluids. This utilizes chiefly the following three effects: (1) By mix- ing the low-temperature cooling gas with the main currents of high- temperature gas, the air-current temperature near the surface of the ^Received March 23, 1971- an First Aerodynamics Department, ***Numbers in the margin ^indicate pagination in the original foreign textf. object Is lowered. (2) The temperature gradient Is reduced on ac- count of the Increased thickness of the boundary layer due to the Increased mass Inside the boundary layer, and also due to variations of the velocity distribution and temperature distribution Inside the boundary layer. (3) Heat absorption occurs because of chemical re- actions or phase changes Inside the boundary layer or on the surface of the substance. Possible heat-protection methods are transpiration cooling, film cooling, local mass Injection cooling, or ablation cooling. They each have their advantages and disadvantages, and It Is difficult to determine which of these methods would be most effective under all possible conditions. Nevertheless, the ablation cooling method con- tains the above-mentioned effects of (2) and (3), and there Is no doubt that It still remains one of the most promising heat-protection modes as long as one Ignores the actual, economical problems such as repairs of the eroded surface of the object. It should be noted that this method has a rather high cooling effect. This method was based on the idea of protecting at all costs the most Important parts of the craft, while allowing, or rather posi- tively encouraging, the melting of the surface substances of the craft by means of external heat. Research concerning this was per- formed at the early stages by Sutton [1] and Roberts [2]. During the process of ablation, generally speaking, a liquid phase boundary layer and a gaseous phase boundary layer are formed on the surface of the object. That is, the solid phase surface of the object melts /2_ producing a liquid-phase boundary layer, over which a gaseous bound- ary layer is formed consisting of a mixture of the evaporated gas and the external gas. The thermal energy in the high-temperature air currents is transmitted to the surface of the object through these boundary layers. Part of it is absorbed into the latent heat of evaporation and the latent heat of melting, while the remainder^ — is absorbed in the form of the thermal capacity of the structural materials. = The substances chiefly used as the surface materials are high polymer plastics, graphite, etc. Other materials are also known In which protection from heat Is attained by bringing about a chemi- cal reaction on the surface or by creating a char layer Inside the object [3, 4, 5, 6, ?]• If the object Itself Is subjected to sublimation, the surface substance will change directly from a solid Into a gas . Therefore, there will be only a gaseous boundary layer on the surface of the object. Nevertheless, this does not Interfere essentially with the validity of the analytical method, which Is based on a two-phase boundary layer Including also a liquid boundary layer. Thus, this paper will deal with sublimating ablation and will clarify Its mechanism. The writer and his colleagues already on a previous occasion used the micro-disturbance theory to analyze sublimating ablation of hemispherical objects, using Teflon materials [5]. The results showed that the calculated values differed considerably from the test values when there was a relatively low stagnation point tem- perature. It would be significant to resort to an electronic com- puter and use a unified method of obtaining a direct, numerical solution, not only in order to check out the results obtained pre- viously, but also to give the theory a more universal validity. From this standpoint, the questions dealt with in this report may be sum- marized in terms of the following three points: (1) Essentially, the ablation phenomenon ought to be determined by the process of chemical reactions of the surface substances in combination with the aerodynamic heat transfer process in the circum- ference of the object. Thus, if the conditions inside the substance and the conditions of the main current are given, the physical vari- ables on the surface of the substance, particularly the ablation speed, the surface temperature, and the heat transfer coefficient, will be determined categorically. (2) The flow becomes dissimilar downstream from the stagnation point of axially symmetrical objects. The; physical variables in the flowwise direction are sought by applying the Hartree-Womersley method, proposed by Hartree and Womersley [8] and further developed by Smith et al. [9] for the dissimilar terms in the boundary layer equation describing this area. The purpose of this method (here- after abbreviated .as| the H-W method) is to reduce a parabolic partial differential equation to a corresponding ordinary one by means of a kind of Implicit scheme .[ The least square error method is applied to the two-point bound- ary value problem of the ordinary differential equation obtained by the H-W method as described above, and a general solution is proposed. A simplified version of this method has been used the Nachthelm and Swigert [10] in cases when there are as many as two unknown boundary conditions. However, in the present case, seven departure values are necessary to carry out numerical integration of the ordinary differential equation. Of these, only one Is given explicitly, and the other six are unknown. These six departure values (unknown boundary conditions) are related by three bo ns train t| conditions. The remaining necessary number of boundary conditions are given at the Itnternal edge|of the boundary layer. Two methods are used in problems of this type: the initial value method, and the quasi- linearization method. In the initial value method [11], the departure values (unknown boundary conditions) are suitably assumed, the numerical values are integrated, and iterative calculations are con- tinued until the conditions of the external edgel of the boundary layer are satisfied. In the quasi-linearization method [12], the original equation is linearized, and a non-homogeneous solution — consisting of a lower-order approximate solution — is sought. Both methods give good results when the number of departure values (unknown boundary conditions) to be assumed is small, but there is a tremen- dous increase in the difficulty when the number becomes larger Rur- thermore the presence of iconstraintf conditions on the wall also complicates the problem. The method used here is a further development of the initial value method. No matter how the number of unknowns or the number of restraining conditions ;are| increased, the method can be applied in principle as long as there is agreement with the fundamental equation. It is also believed that the method can be applied to other problems as well, as long as they are boundary layer type problems (with asymptotic behavior at the upper limits of the domain) , Let us assume a semi-stationary ablation, in which the surface sweepback velocity, and thence the ablation velocity, is almost un- affected by the time. (For a detailed discussion, refer to [5]-) The thermal properties and physical quantities of most ablation sub- stances are not clear [7]. Teflon (polytetraf luoroethylene) is one of those substances which are subject to sublimating ablation (it has been observed experimentally that there is probably a melting process, but since the temperature range of this process is extremely small], it is assumed that there is only a sublimation process) and for which the physical quantities have been clarified. Thus, Teflon is used as the model in these calculations. In the method described here, as long as the state of the main current and the state inside the object are given, it is possible to carry out automatically the above-mentioned processes (1), (2), and (3), and the ablation physical quantities can be calculated in a unified manner as far as any point downstream from the stagnation point of a blunt-headed object. Such methods of direct solution usually rely on the memory capacity of a computer, but in the H-W method the upstream solutions are discarded as soon as they are used. This makes operations less difficult, and this is one of the advan- tages of the method. 2. Notation a: velocity of sound C: Chapman- iRunf sln| number 13- Cp: specific heat at constant pressure Cv: specific heat at constant volume D: two-component diffusion coefficient d: boundary layer thickness E: square error In two-point boundary value problems H: total enthalpy function H „„: effective ablation heat ef f K: concentration function L: latent heat of sublimation M: molecular weight M : Mach number of main stream ! cx> ' m: ablation velocity Pr : Prandtl number p: pressure . - - - q: heat transfer rate R: gas constant Rb : radius of curvature of object Re: Reynolds number r^ : radius of curvature of cylinder Sc : Schmidt number (s, n ) : transformed coordinate system T: temperature t: time (Uj V): velocity components In spatially fixed system of coordinates (u, v) : velocity components In system of coordinates fixed to object V, : surface sweepback velocity due to ablation (X.j.y): spatially fixed system of coordinates ^ ; (xj y): system of coordinates fixed to object 'i X.: boundary conditions z: transformed concentration function specific heat ratio of air error in two-point boundary value problem criterion of square error heat transfer coefficient cooling effect due to ablation coefficient of viscosity density transformed i g^tream functionl stream function] transformed enthalpy function Suffixes: a b e , g n no s sh t w 1 ( )' ( ~) _ K ( ) state with ablation state inside the object state at external edgejof boundary layer state of the gas phase n^^ state state' with no ablation state of solid phase state immediately after shock waves state at stagnation point of main strea: state on wall of object state at stagnation point of object state of ablation gas state of external air state of main stream] differential dimensional quantity quantity for non-dimensionalization Basic Equations Let us consider semi-stationary ablation, in which the surface "sweepbackl velocity is more or less uniform and does not depend on the time. It is assumed that all variations caused by ablation are J7 confined within the boundary layer and have)no effects on the exter- nal current. There are no chemical reactions Inside the boundary layer, and the boundary layer consists of a two-component mixed gas (the gas caused by ablation and the external air). The thickness of the boundary layer is assumed to be sufflclentr y small in comparison! with the radius of curvature of the object. In a flow field such as that shown in Figure 1, the compressible laminar boundary layer equation for a hemispherical object (Appendix A) can be written as below after transformation of the independent variables by means of the Less- Dorodnitsyn transformation: ;; Shock waves Shock wave, layer] 5= I pelitu.rMx, Ji=p,u,r,{2s) * I £-dy Jo J op, '(C5S")'+/3(-e.'-(5>'+l)2}+Vi"(5> + ,)=i?i, + Wu,' |C(1-A)5>"(?>' + 1)}' =i?2, ^&*'^'+*'('*''"'')=^'- (3 1) (3 2) (3 .3) (3 4) Te{X),Ue(X) 'Boundary layer . M., n^ Figure 1. System of coor- dinates and flow pattern /4 Here, ( ) ' is the differential re- lated to n, and a=.l-CNCM^^^ (3.5) C, Pr, and Sc refer to the Chapman-Rube sin number, the Prandtl num- ber, and the Schmidt number (Appendix B) ; they are each functions of (|) , ^3 and z. Here, the values of cj)', ^, and z are the following f=u/u.-l, V^H/H.-l z=K (3.6) The stream function fi is defined as follows:] 3» _„„, 9a fl=V 2s {^is,,,)+,,^ (3.7) Therefore, cj) is expressed as follows: =i=r' ^ / (3.8) If we use an equation of state for a perfect gas, the following will apply because of the condition that there are no pressure vari- ations in the y direction in the boundary layer: f-={i- ^^-^>^}4r.| (3.9) - Due to the equation for the total enthalpy (A.l6), the temperature will be T _ !r+l-|«.W(^'+l)' l , In the basic equations (3-2) - (3.^), the dissimilar terms R-, , Rp, and R^ on the right side will each be: (3.11) j?,=2.{(^'+ i)^i^-r'ii|±£)} J ■( 3 : 1 2 ) (3.13) The dissimilar terms will be expressed by these items and by the differential term of s. The treatment of this will present one of the problems. That is, because of the presence of dissimilar terms R-, , Rp, and R^, it will be impossible to express the flow field by means ;of o ne s imilarity! parameter n , and a separate solution must be sought for each station (that is, for the various values of s). 4. Boundary Conditions Since the basic equations (3.2) - (3.4) are, respectively, three-step, two-step, and two-step differential 'equations with^Te'- spect to ()) , ^i*, and z, altogether seven boundary conditions are neces- sary. However, only four conditions are explicitly given: the con- dition that there is no slipping on the wall (u, = 0), and the W )' condition that there is connection with the external flow at the end of the boundary layer (u ^ u^, H •> H^, z -> ; y ^ y^)- When these are written in transformed form, they are as follows: ,= ; 5i'=-i. I (4.1) ^ = oo ; 5>'=0, 1^ = 0, 2=0 I (4.2) However, the following thre]e Icons traint i conditions occur on account of the physical conditions on the surface of the object where ablation takes place. (i) Relationship between, wall temperature T^ and mass loss rate m When an ablation substance^ is_ considered, __lts_mass-_loss__. rate \m is a characteristic function determined by the wall tempera- ture T and the pressure p. It can be written as follows: w m = FCTu,.J)\ (4.3) However, various test results and data [13] have shown that in the low-temperature region the effects of pressure are so small as to be almost negligible. Rashis and Hopko [l4] have obtained this relationship for Teflon as shown in Figure 2. When this is approxi- 75. mated in a quadratic equation, it becomes: m = — ptVi= (pv)w (4.4) One can thus obtain an approximate equation to take the place of equation (4.3) • On the other hand, from the aerodynamic viewpoint, if the state of the main ■streamj is given, the relat ionship between the wall tem- perature T and the mass loss m can be described by the]- .. following equation: ^=G(T,o, 0)1 (4.5) 10 7.0 S56.0 E V5.0 ^4.0 3.0 2.0 1.0 0.0 « * 5.0 Si 4.0 E (J ,3.0 X IS 500 550 600 650 700 7501 ru<'c) 2.0 1.0 /?6= 0.75cm i? t= 2.00cm 3.00 Figure 2.. Chemical rela- tionship between the wall t emp e ra ture an d the ma s s loss rate for Teflonl [TiT] . 500 550 600 650 700 Tw{'C) 750 Figure 3. Aerodynamic rela- tionship between the wall temperature and the ablat-lon rate for Teflqnj Here, Q Is a quantity which Is determined according to the Mach number of the main ,s treaml , the state of the. stagnatJ:On point of the T andm w are ' unknowns main stream), the state of the object, etc In Equation (4.4) and Equation (4.5), and It Is anticipated that the two unknown numbers will be determined uniquely with respect to various ambient conditions by converting these two equations Into a simultaneous equation. In actual fact, the relationship in Equa- tion (4.5) was obtained as in Figure 3 according to [5]- By over- lapping Figure 2 and Figure 3 diagrammatlcally , it is possible to obtain the solution with T and m at the respective intersection. w In this report, we proceed one step further and attempt to—at- — tain a comprehensive solution, in which Equation (4.4) is also in- cluded as part of the boundary conditions. When Equation (4.4) C 11 is written In transformed fashion. It assumes the following appearance: X {Ai 4- AzT.T ,0 + Ja (J I r,o)2J (4.6) (11) Mass equilibrium on the surface of the object K. , the concentration of the various component gases on the iw surface of the object, cannot be designated without relation to the -Wall velocity v . When v is determined, k. will also be deter- f| WW IW mined. Generally speaking, v^ , the velocity of the 1 gas in the y D dK,l direction, is the sum of the diffusion velocity "a", ITjand the con- vection velocity v. Therefore, pK,v, = -pD^J^+pKiV dy (4.7) Only the gas components produced by sublimation will be blown awayj from the surface' of the object. Therefore, one can posit the physi- cal condition that there is no movement of the external air Inwards or outwards along this surface. Under this condition, the v^^ with respect to the external air (suffix 2) will be 0. If we use the fact' that Kp = 1 - K, we obtain the following: Pvo H-f)- H-p. V. d-K", )=0 (4.8) Equation (4.8) can be rewritten as follows in transformed fashion: m. r^=^-2.(^), (4.9) (1). Ata constant stagnation point temperature, if there is an increase of T there will be a decline of the temperature gradient inside the w boundary layer. Consequently, the thermal input towards the surface of the object will decline, and there will be a smaller ablation velocity. 12 (ill) Energy equilibrium on the surface of the object On the surface of the object where ablation Is taklijg place. part of the heat quantity V dyJz which Is about to be transferred through the boundary layer to the object will be absorbed by the latent heat of ablation, and the remainder of It will be transferred Inside the object. Thus, If we take Into account the conduction of heat In the x direction Inside the object, the energy equilibrium on the surface of the object can be written as follows: ('41. =^^>~>-+-{(^);.. dy /». (4.10) The second Item on the right side Is obtained as the solution of the heat conduction equation for the Interior of the object: _,3'T fT\ dT_ _:._ n^ (4.11) Nevertheless, since Equation (4.11) is an elliptic partial differen- tial equation. It would be Impossible to combine It In a simultaneous equation with a parabolic boundary • layer equation and apply the H-W method. Furthermore, with the range of stagnation point temperatures being handled here, heat diffusion In the x direction Is negligible In comparison with that In the y direction [5]- Therefore, Equation (4.11) can be written as follows: dT dy (4.12) Here, m = - p^ v^ = ( p v)^ /6 Let us solve Equation (4.12) under the following boundary conditions : 7=0 ; f = T„ -1 / y=-yr.T^Tjj (4.13) In dimensionless form, the following solution will be obtained: r=c, exp(^^)+J (4.14) 13 Here , c,= (r.-r.)/{exp(— -|L^.)-i|^ Cz=T„—ci, R VS. (4.15) Therefore, the following will apply dT\ _ P (3T\ P l\-l. (4f)., .=P' /""•''•' (4f)., . (4.16) By substituting Equation (4.16) in Equation (4.11), we will obtain the following equation for energy equilibrium in transformed fashion: 1 — aZm ' PalCtoUtTeTt Here {-^ gg-/+W^-+M (4.17) E\w /=Z/C,2'5",. (4.18) The wall temperatures T appearing in the various equations can be expressed as follows, in accordance with Equation (3.10): r. = (!P-. + l)/(l-a2.) (4.19) They are related to ^ and z . ^ WW 5. Calculating Methods 5.1. Hartree-Womersley Method The H-W method, which was described in the Introduction, trans- forms into the difference system the dissimilar items R, , R„ and R on the right side of the basic equation, which is a parabolic partial differential equation, as well as the differential terms] having the 14 form of dF/ds which are contained In the boundary conditions. The equation Is thus reduced to an Implicit ordinary differential equa- tion. The two-point differences of equal Intervals used here were: OF _ F-F.-i l (5-1) ds As / However, the precision will Improve further If three-point differences with variable Intervals are used. Nevertheless, the conve^lon velo- city distribution on the surface Is smooth In ablation problems. Therefore, there Is no necessity for discussing the difference sys- tem with reference to the discontinuous convection velocity which was dealt,, with In [15], and It will be sufficient to use two-point differences . According to the calculating diagram In Figure 4, the solution for the n point (the suffix n Is omitted In order to make the equa- tions easier to read) will Include the solutions for points upstream u n-2 n—1 from It. That Is, at the nth point, one will use F„_-| i the known solution for the points upstream, and the ordinary differential equa- tion for n will be solved in the n direction. The dissimilar terms on the right side of the basic equation will be written as follows : Solu-I tlon Fn-ii. 1)) Solu- tion Fn-iiv) Solu-] tlon I ■^S Figure 4. Calculation dia- gram for the H-W method i?i=-|^{(i' + l)(*'-0'.-l)-«>"(!*-?>.-l)} As mentioned above, 4 -,,(!)' -,5*1' -, , z -, are known solutions ' ^n-1' ^ n-1' n-1' n-1 (5.2) 'The stagnation point of a blunt object becomes the point of departure for calculations, and there are no solutions for points upstream from it. However, with respect to the finite difference amplitude at this point. As, it can be stated that 2s/As = 0,. and therefore one can obtain a solution corresponding to cases where Rj = R2 = R3 = (similar solution) . 15 5.2. Two-Point Boundary Value Problem The ordinary differential equation obtained above does not have the necessary number of boundary conditions on the wall, and It therefore becomes ja two-point boundary value problem of the type described in the Introduction. Let us posit the following in order to give the problem a uni- versal validity. "In a boundary layer type ordinary differential equation with /7 n steps (asymptotic behavior), m (m < n) boundary conditions are given at the .lower limit of the domain, and (n - m) of the conditions are unknown. These (n - m) unknown boundary conditions are linked by means of I {I < n - m) constraint] conditions. At the upper limit of the domain, (n - m - I) boundary conditions are given. We seek a solution satisfying these conditions." Let us express the boundary conditions In terms of Xj (j = 1, 2, ..., n) . First of all, the conditions which are explicitly given are X-,, Xp, .... X (their number is m) . In the initial value method, (n - m - I) boundary conditions, excluding this and the I constralntj ' conditions, must be supposed as the departure values for numerical integration. These boundary conditions which are to be supposed as expressed as X^^^-j^, \+2' "' "^n-l ^^^^^^ number is] (n - m - Z)]. If these (n - m - I) conditions are suppos ed, the remaining I conditions will be determined from the ^c ons traint conditions, and therefore the solution insidje the domain will be governed by the Independently assumed! X , -, , X ^„, ...., X ,. If the solution is expressed as S., - - — I m+1 m+il n-i J it can be written as follows : Sj =Si (Xm-ny x„+2, , x.-d \ (-5.3) (y=i. 2, , «)1 16 In the present problem, n = 7j m = 1, and I which are given explicitly are: and the conditions which are to be assumed] are ^2=5>"^, Xi=W,, Xi=z,\ 3; the conditions (5.5) The conditions which are determined via the ;co-nstralnt| condi- tions are: Jf6=/. =X6(Jt2, ^3, X4) = ( 4 • 9 ) ^ (5.6)/ The solutions inside the boundary layer are: Si=55=Si (X2, X3, X^) S2=55' = S2 (X2, Xz, Xt) S3^p" = S3(X2, X3, X'4) St^'F = St (X2, X3. Xt) S6=r=S5(X2, X3, Xt) S6=2=S6 (ATa, A'3, Xt) S^^^' = St (_X2, X3, Xi) (5.7) As was'' mentioned before, the number of conditions necessary at the end of the boundary layer is (n - m - Z). In the present case, this number would be 3. They are given by equation (4.2). However, as was mentioned in the Introduction, in these boundary layer type equations, all the physical quantities have an asymptotic behavior at n ■> «>. Therefore, the following also ought to apply supple- mentarlly: l? = oo ; ?>" = 0, if' = 0, i'=0 (5.8) If the conditions in Equation (5.8) are not present, there may be solutions in which ^a certain quantity will overshoot along the way (Figure 5), so that the solution is no longer an asymptotic on e. In order to avoid this phenomenon, the supplementary conditions in Equation (5.8) become necessary as full conditions. 17 Let us use suffixes I and II to indicate the first approxima- tion and second approximation of solutions at the external edge of| the boundary layer (the upper limit of the domain) by Iterative calc ulations . ^jln this case, the modified quantities of the de- parture values (the supposed un- known boundary conditions and the boundary conditions determined by the ' constraint! conditions) and the modified quantities of the solutions will be: Figure 5' Explanation of solutions, with overshoot in the center ^X,= (,XOii-(iXi)i, (.-=2, 3, 4) ^S/=(Sy )ii -(Sy ) I , (; =2,3, .7) Furthermore, according to the yairlatlonal principle ASy=2]^AXf (;=2, 3, , 7) (5.9) (5.10) (5.11) Here, the operator D DX. 1 is defined as follows : ^ ydXk d DXr'dXi'^LdXi dXk According to Equation (5.10) and Equation (5.11) » 4 '(Sy)lI=(Sy+2]^^X)l (/=2. 3. 7) (5.12) (5.13) Since Equation (^.2) and Equation (5.8) are the desirable conditions at the external edge of the boundary layer of (S.)jj, If <S . is the error between (S.)jj and them, it will be: '-Sj =Sy +J]^AXi U = 2. 3, , 7) (5.1^) The suffix II was omit'ted here in order to give universality to the approximations. It is clear that the necessary and sufficient condi- tion for convergence of the solution is to give 6j a value of 0. For this purpose, let us take a square error of /8 18 L^' (5.15) and assume that the convergent solution will be when this assumes a value of 0. In actual practice, let us take a small positive number and adopt the following as the criterion for convergence: £ < £ (5.16) The conditions for the least square error are obtained from the following: dE d (AXi) = 0=2, 3, 4) In the present problem. dE 1 = 2 1=2 4 g-^=0J:Oe„ + ^e.<AX<=0. (5.17) (5.18) Here y=2 7 »-i:(S>' . (1 = 2, 3, 4, *=2, 3, 4) (5.19) The following is obtained from the above simultaneous equation: 19 AA'2= - AX3= - AA'4= - DEM = ©21 023 024 031 033 034 041 043 044 ©22 021 024 032 031 034 042 041 044 022 023 021 032 033 031 042 043 041 011 012 013 021. 022 023 031 032 033 /DEM , /DEM /DEM , (5.20) In one trial, a modifier term of AX^ is added to X^, giving X^ + AX^ , and we move on to the next trial. However, in order to obtain AX. in the manner of Equation (5.20) it Is necessary to find, by some method or other, the value in Equa- tion (5.19) taking the form of DS./DX. . This is necessary at all the stations in the boundary layer. Differentiating the basic equation by X. (1 =2, 3, ..., 7), we obtain the disturbance differential equation : / s ^+^A{-^-('''-^»^&^'> + !»' „S^_dR\ 3 f C ,\.3z',., s ,,, # _dK3\ («=2, 3, •. 7) =0 (i^j) («, y=2, 3, •. 7) (5.21) (5.22) (5.23) Solving this equation under the following boundary conditions : (5.24) let us substitute this solution, together with the values of 8X^/8X^ (1=2, 3, 4, k=5, 6,7) (known on account of the constraint 20 conditions between the boundary conditions). In Equation (5-12) , the defining equation of p-^:— . The S. In Equation (5.2^1) Is a function determined as In Equation (5-7). All of the boundary conditions (departure values )' for the disturbance differential equation obtained here, are given as In Equation (5.24). The solution Is obtained by making a simultaneous equation of the total of 21 differential equa- tions. Including the 3(n - m) disturbance differential equations (there are l8 of them In the present problem) and the three original equations. The square error Is calculated at each of the calculating points (at each of the scale Intervals of the numerical Integration In the n direction). Nachthelm and Swlgert [10], established In ad- vance a suitable n a. (point for stopping calculations) and Investl- StOp -^ Ir-Jr- o gated the effects of Its value on the solution. In this report, we posit the following criterion: |Sy(,)|^a,y| ( 0) . = constant) (5.25) J Solutions departing from these conditions are regarded as divergent solutions, and the calculations are stopped at that point. This point is called n = n ^ • Then we proceed to the next trial. The ^ stop end point of the calculating region given in advance is also called n = n . . In this way, the convergent solution is one which matches ' 'stop '^ ' ° the conditions of Equation (5.I6) while also fulfilling the condi- tions of Equation (5.25). Thus, the n = n satisfying the conditions of Equation (5. 16) will correspond to the external edge| of the bound- /g ary layer, and it is possible to find also the thickness of the boundary layer, which was one of the unknown numbers. The flow chart for the calculations is as shown in Figure 6. 21 x = x+/ix C 111 ,'") s = s +Js :« ;! ( »• = 2 , 3 , 4 ) >7=0 3?= J7+ JJ? (Runge-Kutta-Gill) El^ _('x. \3^ J-X-f ( 1=2,3,4) V=ve ^:3Z7 Figure 6 1 — beginning; 2 Flow chart for the H-W method and the method of solution for two point boundary value problems - read In conditions of main .■jstreamr and Inside the assume Initial values of X 3 — stagnation point; - Xs, Xs, and X? are determined from lconstralnt|, 2.1 J X 3 5 object ; X:^; 5 — Xs, Xs, and X? are determined from lconstralnt|, conditions 6 —original differential equation and disturbance different lal_ equa- tion are solved as a simultaneous equation; 7 — - E Is calculated; 8 — E Is written out; 9 — end; 10 — does the solution converge?; 11 — (downstream region); 12 — obtained solution is included in Ri , R2 , R 13 — solution is written out; l4 — AXj_ (1 = 2, 3, 4) is calculated." 3 5' 22 Numerical Calculation Results and Considerations of Them The Teflon (polytetraf luoroethylene) used in the numerical cal- culations is a fluoride polymer, but it decomposes into a monomer upon pyrolysis. According to Madorsky [l6], more than 90^ of it is CpF. . Therefore, let us here regard C2Fij as the ablation gas and assume that the boundary layer consists of a two-component mixed gas (air and CpFh gas). The physical values of these substances are as follows : C2F4 gas: ^^^^00 j Cpi=0.32 cal. gr-i. ."K-i/ ffi=5. OOA r s j2,(2.2)=o.90 / (See Appendix B) /lO Air: A/2=29 Cm=0. 28 cal . gr-i . °K-i 72=3. 62 A .O2«-2)=0.76 (See Appendix B) Solids Teflon : c>.=o.22 cai . gr-i. "k-* pi =2. 19 gr .cm-a Kt = 6. OOX IO-*''cal.cin->.sec-». °K-> L=35kcal • mol-i The following are assumed as the conditions of the main stream : A/oo=5.74, r,=800°C~1500°C, r.=20''C, ^( = 0. 5atm~2atm, ^»=0. 75cm~3. QOcni In Figures 7(a) and 7(ti) is shown the state of convergence of the departure values of a two-point boundary value problem when— the— calculations were performed according to the method described in the preceding chapter. The calculations were carried out on the HITAC- 5020 computer at the National Aerospace Laboratory. Figure 7(a) shows the state of convergence between (j)'' , W m w ^w' the values /ll ^ 23 -1.0 ¥m -0.8 -0.5 z'w --0.4 |l-0.6 ii 1 ^0.4 1 — r 1 — I — I — r , -0.2 -0.1 1.0 'w, Zw : Departure value s j ^^'=0.0, 7'»'=0.90,arH-=0.'01 :rlals E 5.41191 1 1.55245 2 2.16962 3 0.11206 4 0.00029 5 0.00008 2.0 p; 1.5 0.5 1.0 0.5 ^y 3 4 5 6^ 7 'No. of trial si Figure ?• State of convergence of solution : (a) convergence of initial value: a; = 0, T, = 1000°C, ^. = 1.00cm, ^, = latm. tO.5 -1.0 Figure ?• State of converg- ence of solution: (b) velocity profile, i*=0, T, = 1060°C, ^ = 1.00cm, i>;=latm.| assigned at random at the beginning (number of trials = 0), and d) , 1" , z' , determined according to them. The square error E in each trial is also shown here. The square error was calculated at n = n • As for the divergent solutions departing from the criteria of Equation (5-25), they are those at the maximum n position (ffT ) S G Op as long as the conditions of Equation (5.25) are satisfied. The following were adopted as the conditions of Equation (5.25): 24 032 = 1.5 (for *') 0)3 = 1.0 (for *") 0)i| = 0.5 (for '^) ^5 = 2.0 (for >!" ) 1^6 = 0.5 (for z) 0)7 = 2.0 (for z') (6.1) The state of convergence of the velocity profile is shown in Figure 7(t) ) . When a small positive number e is determined, as in Equation (5.16), it becomes a parameter for determining the boundary layer thickness d. Therefore, attention is focused on how its value is selected. The error will be 0.01 if the external end of the boundary layer is the point where the physical quantities in the boundary layer amount to 99^ of the main jstreaml . . Therefore, a value of '-B^o.oPx 6=0.0006) -4 ■'is desirable. Thus, a value of e = 10 will be quite sufficient. In actual fact, according to the calculation results, the relation- _ _ _4 ship between £ and d (normalized at d when e = 10 ) will be as shown in Figure 8. 'It is only natural that if a large value of|£ is adopt ed,'ithe boundary layer will be evaluated with a reduced thickness. Even if the value is less than e = 10 , it is known that there will be a considerable reduction. In view of the above considerations , -B^«=io-*) was used as the criterion for convergence. The next question which arises is whether a uniform value can be obtained regardless of the departure values (the uniqueness of the solution determined by the initial values ) . The departure values /12 were varied (cj)"^ = 0.0, 0.25, O.5O, 0.75, 1-00; T^ = 0.7, 0.9; z^ = 0.01, 0.21, 0.41) in order to study the numerical uniqueness. As a result of these calculations, it was confirmed that in all cases there was a convergence to a uniform solution. Typical examples of this are shown in Figure 9. Uniform solutions are obtained even though the number of trials differs. It Is possible to select these departure values at random, but if the values are selected in the 25 1.3 1.2 1.1 1.0 0.9 0.8- 0.7- 0.6- 1 1 — I I I I i i | 1 1 — I I I M i | I I I I I I m r I — I I I 1 1 d/d 6=10"* 0.5 10-* III I ' ' 1 1 -I I I Mill 10- 10-" 10" 10- Figure 8. Effects of the convergence criterion on the boundary layer thickness x=0, 7, = 1000°C, ^»=1.0pcm, ^=latm.| vicinity of the anticipated values, naturally this will contribute to speeding up the calculations . The calculations performed here were for cases when there were three unknown boundary conditions (departure values) and three constralntjcondltlons . However, It Is known that this method Is valid when a similar method Is applied In cases with one or two departure values [15] (the work Is much easier ,. since there are no constra int ! conditions). Therefore, the general argument In §5.2 Is believed to apply also In cases when there are four or more unknown boundary con- ditions. (The task of numerical confirmation Is left for future study. ) As was mentioned In §5.2, It Is desirable to consider the sup- |, plementary condition equation (5.8), as this will give more perfect :, results as the necessary and sufficient conditions. However, even If ., this is omitted, there ought not to be any deficiencies as far as the , 26 necessary conditions are concerned. Assuming that the former is Case 1, and the latter is Case 2, calculations were per- formed under the follow- ing conditions : 1273° K, R^ = 1.00 cm. \ = p, = 1.0 atm, and the de- parture values of (j)" ^ w 0.0 = 0.01, T^ = 0.90 Departure values 0.40 Zic 0.30 0.20 • O^«'=0.00 A rH.=0.70 □ Z«r=0.01 ATh'=0.90 I Gf 2^=0.41 0.10 "W ' w When both cases were com- pared, the results shown in Table 1 (convergence values of the departure values) were obtained. There were no very great differences and, as was only natural, the calcu- '> lating time was less in Case 2. The lack of dif- ferences in the converg- ence values is based on the nature of the solu- tion, as described in §5.2. That is, there are no overshoot solutions (Figure 5)} and all the physical quantities have an asymptotic tendency' m^ving_ towards 0. / All the calculations below were performed according to Case 1, taking Into consideration the necessary and sufficient conditions. 4 5 6 No. of_t rials] 1 Figure 9. jf=0, Moo=5.74, r,= 1573°K, 5,=0.75cm, p, = {atm. Uniqueness of the initial values /1 3 A profile of the various physical quantities inside the boundary layer obtained in this manner is shown in Figure 10.. If cj) and T, ■^ WW 27 TABLE 1. COMPARISON OF CALCULATION RESULTS BY DIFFERENCES IN_THE BOUNDARY CONDITIONS AT THE EXTERNAL EDGE) OF THE BOUNDARY LAYER [Case I: Including the supplementary condi- tions of Equation (5.8)'. Case 2: Not in- cluding the supplementary conditions of Equation (5-8)] Case 1 Case 2. Number of trials 6 il 5 5». ■ -0.02873 * -0. 02895 ', <>'" 0.79907 ' 0. 79922 T. 0.71379 0.71418 z. 0. 05375 0. 05381 z'. -0.02719 -0. 02749 E 0. 00008 0. 00007 V- 3.7 132 s ec| 3. 7 104 s ec| Calculat (using ing tin HITAC le, . 1 5020S)1 are known, it will be possible to obtain the ablation rate m a nd[ the surface temperature T from the following: w * p.fi.u.r.^ Y' , /, ,2s\ 2s ^ -\ (6.2) (6.3) r.=r,r. I Figures|ll(a) and (b) show the effects of the stagnation point tempera- ture on the ablation velocity and the surface temperature, plotted with the radius of curvature of the object as the parameter. When a comparison is made between cases when the thickness of the object was made, infinitely great /(y^; = °°) — ^^ thaf'ls, "wTien the object was a \ solid hemisphere — and cases when the thickness of the object was J flnife "(1/2" of the radius of curvature of the object; y, = 0.5 R, )|^-' •' tha-t--is, when the object was a hollow hemisphere with a thickness . - half the object's radius of curvature — it is found that. in the " • ^ latter case the ablation velocity and the surface temperature are rather low, while the heat transfer rate, shown in Figure 11 (d), increases on the contrary. This is due to the following fact. 28 1.0 T i' 0.9 0.8 -0.5 — ^ 0.7 0.6 .0.5 -1.0 ■■■ 1 -0 1 1 1 I \ ■* -^ / /^ V \ s^ J ■ \ / f -0.1- 1 / I \* / r"" ^_ \ / / ^ \X \ / / \ 1 1 \l fl \ Ik ^ -0.2 - \ / /' V y ' - -\i.o//' ^ <., , v / - yi -0.3- /A // \ / _z_ //' \ /'^^ //' ^i . III \ « "' / \ -0.4 - 1 \ 1-2.0 , L... ^"^^sL , -0.5 0.10 z 0.05 Since a uniform tempera- ture T, was assumed at y = yi_ J in the latter case it is necessary to maintain a lower tempera- ture on the surface of the object and to suppress the ablation effect. The thermal input increases for this reason. In order to indicate the effects of the object thickness, calculations were made of the ablation velocity at the stagna- tion point of the object and the surface heat transfer rate at T, = 1000° C. The results are shown in Figures|12 (a) and (b). These results indicate that, when the object becomes thinner at the same thermal input, there is an increase in the heat transfer rate because of the necessity to suppress the ablation velocity; this results in a deterioration of the ablation cooling effect. As the object thickness y, becomes smaller, the deterioration of the cooling effect proceeds rapidly ahead. Although there is little possibility of using an object which is excessively thin, it is believed that it would be possible to put into actual'ase 5 6 V Figure 10., Profile of the physical quantities in the boundary layer 'x=0, ■T,=1500°C, ■S.=1.00cm, ^=iatm.l The difference between the surface ^b a hemisphere of R = 1.00 cm. heat transfer rate of an object of y, = 0.25 Rj^ '^ 0.5 R^ and that of 717 29 xlO"' '(«Mm-''K<!-') -1 — I — I — r -1 1 1 1 — JT- i 2 - =2.00cm .J I I I I I ' ' ' 500 1000 ISOO 7.CC) Figure 11(a). Ablation rat e x=0, Mco=5.74, ^(=latm. 680 Twi'C) 660 620 580 T — I — I — I — I — I — rrr _1 I I I I u 500 i 1500 : Figure 11(b). Surface temperature *=0, il/oo=5.74, ><=latin.l this method (y^ = «') ; this method (y, =0.5 R,); [5] (y^ = -) 0.15 Z«.- 0.10 0.05 I I I I I 1 1 — I 1 — I r y»=0.5;?» S»= 0.75cm _1_ I I I I I I 1- 500 1000 ■r.co 1500 5.0 4.0 - 3.0 - 2.0 1.0 1 1 1 r 1 I .nfM..,,, |7. rn 1 1 1 ■ »» = 0.5R, y?.=0.75cm / 1.00 ^ f /" ^ 7 / ■■\ — / ; '^ 2.00cm 1/ /'' y "^\^ 3.00 500 1000 r,ca 1500 Figure 11(c) Concentration of Figure 11(d). Surface heat trans- ablation gas on surface of object fer rate when there Is a blation ir=0, Moo=5.74, /i7=latin.) «=0, Moo=5.74, A=latm.] 30 1.1 1 1 1 1 1 1 1 1 -1 1 1 1 1 1 - 8 ,'iio - ^*=oo __^-.— . O <0.9 ■i 0.8 "t^ V\\/?6= 1.00cm W 2.00 \ 3.00 0.7 - ] - 0.6 - - 0.5 n 1 'i 1 " 1 1 1 1 1 1 l.D 2.0 3.0 4.0 _ _ 5.0 VblRb Figure 12 (a). Effects of the thickness of, the ■ object on the ablation rate • '■'■'- ^ 1.06 r Qo/(Qo) 1.05 1.04 1.03 1.02 1.01 1.00 *=0, Afoo=5.74, r,= 1273"'K, />,=latni., -'-^ '■■'■■'■ '■■'■■T-r-1" r 1 — 1 ' - - -1 ■ ^(,= 1.00cm \ y 2.00 - I // 3.00 J ^^ • * - ' ' ^^ y;=oo \ • 1 1 1 1 -« J 1 1.0 2.0 3.0 4.0 _ ,5 5.0 Figure 12(b). Effects of the thickness of the object on the surface heat trans- fer rate x=0, M=o=5.74, f;=1273''K, ^=latinj 31 0.4 0.3 0.2 0.1 T — I — I — I — I — I — I — I — I — r gt=3.00cni T^ 2.00 1.00 7^ 0.75 i?»=0.5/?» J \ I I I I I L_l I I ' 500 1000 1500 TA°C) Figure 13Ja)_; Heat protection effectl,of ablationi i=0, Moo=5. 74, /'< = latm. 1000 (calgf-') f.CC) Figure 13(b). Cooling effect of ablation jr=0, M» = 5.74, /!r=latm, J,= oo| an object of y, = °° (a solid hemisphere) is only about 6^ at most. This is attributable to the good thermal insulation properties of the ablation substance Itself. Figure 11(c) gives the concentration of the ablation gas on the object surface, which was calculated at the same time. It is clear that there is the same tendency as with the ablation velocity. Figuresil3(a) and (b) show the cooling effect of ablation. First, Figure 13(a) shows the ratio between the heat 'quantity ^ ~ ^h which is protected by ablation and the heat quantity q when the-re — - is no ablation. It is clear that it increases together with the increase of the stagnation point temperature of the main fistreamjj]f| and ' that the cooling effect improves at higher temperatures. The fact that the cooling effect is better in a solid hemisphere than in a 32 hollow one points towards a reasonable result, namely, that a larger ablation effect gives a better cooling effect. If we define the cooling effect of ablation as : A=-^^ Hw fft9 (6.4) the value of this will display results of the same type as those above, as is clear from the broken lines in Figure 13(b). When cal- culations are made of the effective heat of ablation, which is often used as an index of the ablation effect [17] 3 H,„=—^ (cal.gr-i) ffto (6.5) this also gives the same results as those for A, as is shown by the solid lines in Figure 13(b). The fact that there is a superior cooling effect when there is a bigger radius of curvature can be explained in terms of the follow- ing circumstances. The boundary layer in this case becomes thicker. and there is a larger heat quantity to be absorbed] inside the boundary layer -i with the^identTcal heat input] ' from the main stream'.] The slight differences from the results reported in the literature [5] are attributable to the manner of evaluating Equation (4.4). If Equation (4.4) is further refined, it is thought that this difference will become smaller. The uniqueness of the solutions obtained from this direct solution method is guaranteed also by the fact that there are no big differences between these solutions. Ablation is accelerated when there is an increase in the stagna- /I 9 tion point pressure of the main current. In Equation (6.2), since ^ "I V is in direct proportion to the - -^ power of p , [refer to Equa- tion (A.23)], the following will apply if the other conditions are uniform: m«i^» «r*«'^.»^[ , Since p.»«^| , in the final analysis the abla'tToh ^_ velocity will be equal to the 1/2 power of the stagnation point pressure. 33 The first problem arising In numerical calculations by the H-W method Is the effect of Ax, the step] amplitude in the x direction. This effect Is shown In Figure l4. The ablation velocity at various points (x = 0.20, 0.40,. 0.60) was calculated for three different step] amplitudes: Ax = 0.20, 0.10, and 0.05. The results were compared with the hypothetical ablation velo- city at Ax ^ 0, and It was found that there was an error of only about 0.3^ at the most at Ax = 0.05- As can be calculated by Equation (3-1). When Ax = 0.05, there will be the following value for a' step 1.1 thilHa 1.0 0.9 - amplitude In the n direction of n = 0.1: 0.8 1 — 1 1 r—" X =0.20 / 0.40 1 /x =0.60 1 1 1 0.05 0.10 0.15 0.20 Ax 0.25 A./(A,)2<M3^^ (6.6) Figure l4 ;ude 7i = 970°C, ^.= 2. 83cm. Effects of stepl amplitude Ax^ln H-W method This also fulfills the stability conditions for the difference solu- tion of a parabolic partial differential equation and Is believed to be suitable with reference to the numerical treatment. Since s = 1 4 •jj- X near the stagnation point, it is true that As = 0.003 even when 2 Ax = 0.10, and the stability conditions are satisfied when As/(An) 1 0.3 < 2* However, the conditions become unstable downstream from this. In view of the above results, the calculations were carried out with Ax = 0.05. In Figures 15(a), (b), and (c), the ratio between the ablation rate ml(x) in the flow direction calculated in this way and that at the ^stagnation point is shown. The test values and the solutions obtained by the microdisturbance theory derived from [5] are also shown for the sake of comparison. The present solutions were seen to have an even better correspondence with the test values. Stable solutions can be obtained only up/ to about x = O.Tj and it is /20 34 1.2 I 1 1 \ 1 r- Mlcrodlstu rba n ce| theory I'd^Tl 0.5 'This metho_dl o Test_ value [5] | J I I I L 0.2 0.4 0.6 0.8 1.0 1.2 X 1.2 1.1 tk/tko 1.0 0.9 - 0.8 - 0.7 - 0.6 0.5 1 1 1 1 1 -^ c -^-5-lJ ■ - Thlsll^\N/\ •method! /^\ \ m=m„cos(fjy^ \ \ o'Test value [5] . 1 1 1 r 1 0.2 0.4 0.6 0.8 1.0 1.2 Figure 15(a) . Vari ation of] of ablation ve'locity (a) iWoo = 5.74, f.==970°C,- 1 Rt=2.S3cm, p,= latm. Figure 15(b). Variation of V of ablation velocity (b)Mc»=5.74, t;=1045°C, S"»=1.40cm, ^, = latm, aj;=0.05 impossible to obtain any solutions downstream from this point on ac- count of the essential validity of the hypothesis and the accumulation of numerical errors . The model used in these calculations is that of a solid hemisphere (y, = <») . The hypothesis of semi- stationary ablation which was used in developing the theory holds that ei'^her the amounts of the varia- tions do not depend on the time or that the time-clependence is ex- tremely small. The ablation velo- city may justifiably be considered a function of x solely. As is 1.2 It I 1.1 1.0 0.6 0.5 T 1 1 r Thls_method| [5] m = TWoCos o.Test value [5]| J -. I I I L 0.2 0.4 0.6 0.8 1.0 1.2 X Figure 15(c). Variation of of ablation velocity (c) Moo=5.74, 7-«=1120''C./?. = 1.15cm, ^, = latm, 6j:=0. 05. 35 shown In Figure l6 , It has been observed in experiments [l8] that, after the elapse of a definite time, the relationship between the surface sweepback distance 7] and the time t becomes more or less linear, and that moc-^ Is either unrelated to the time or has only (2) an extremely small dependence. ^- In the case of a perfectly steady ablation r|l9], the following would apply: m(i-) = m.co9(-|^)j (5^7) /21 The values are shown graphically for various different cases. The slight differences are due. It Is believed, to the dependence on time. Within the range of lower stagnation point temperatures In the main .stream] , the thermal Input Itself Is found to be low In the experiments, while there Is little heating at all points except those In the vicinity of the stagnation point In the object Itself. There- fore, the distribution of the surf ace__s_we_epb_ackJveloclty becomes smaller as one moves downstream from the stagnation point, and It becomes more and more difficult for the seml-statlonary hypothesis to remain valid. To sum up the results of the calculations, the ablation velo- city at X = 0.6 can be expressed as In Figure 17 with the radius of curvature of the object as the parameter. The distribution In the X direction of the heat transfer rate and the surface temperature are also calculated In a similar manner. The question at issue here is the conduction of heat in the x direction Inside the object. Within the range of main current stag- I natljDn point temperatures handled here, this has a value which is negligibly small in comparison with that in the y direction [5], but within the higher temperature ranges it becomes doubtful whether this is valid or not.. Popper, Toong, and Sutton [20] have recognized ^^^The fact that I"! Increases in the negative direction at first in Figure 16 indicates that there is elongation due to the thermal expansion of the model. 36 6.0 ^ 5.0 4.0 i 3.0 2.0 I.O 0.0 -1.0 ■ — r Mode \ y ..^y<X\ V. . y. <y V-^y^/^ \ :/ y ^ — -^ 10 1 1 150 200 250 r(sec) Figure l6. Measured values of surface sweepback distance for Teflon [18] ;Ma.=5.74, f, = latm T,= 1200''C, ^=1.21cm.| hlrria 1 1,1 1 1 1.1 ^6= 3.00cm 1 n / 2.00 - III 0.75 0.9 ^^^====^r^ 0.8 ^^^^^^^ ^^ 0.7 ^ ^> — N X 0.6 II 1 1 1 1000 1200 1400 T.CC) Figure 17- Ablation velocity at X = 0.6 A/<»=5. 74, ^ = latm. the results obtained with the mlcrodlsturbance theory of the writer and colleagues within the low temperature range and have obtained good results with tests at higher temperatures using beryllium and graphite as the materials. Their tests Included the three-dimensional heat conduction effect. 7. Conclusion The direct solution method for sublimating ablation Is one of the most effective methods of thermal protection for objects re- entering the atmosphere at hypersonic speeds . This method was studied In a unified manner from the following three standpoints. \ (1) Essentially, the ablation phenomenon ought to be deter- mined by a combination of the chemical reaction process of the sur- face substance and the aerodynamic heat transfer process in the ; 37 circumference of the object. The numerical findings obtained here indicate that, if the conditions inside the object and the conditions of the main streamj are given, it will be possible to solve the boundary layer equation by adopting the chemical properties of the surface substance as the boundary conditions on the surface of the object and also by using the heat conduction inside the object as the boundary conditions in the form of the energy balance on the surface of the object. By solving the equation in this manner, it Is possible to obtain categorically the physical quantities such as- the ablation velocity or the surface temperature. (2) The Hartree-Womersley method was applied. In this method, the differential term of one variable in a parabolic partial differ- ential equation is replaced by a difference term, and the equation is reduced to an ordinary differential equation. The solution for points up to X = 0.7 downstream from the stagnation point was ob- tained as in (1) above. Since the distribution of the ablation velo- city in the x direction is smooth, the stability conditions are not very rigid when the differences are taken. ^ -- •- (3) Since the number of boundary conditions on the wall' does not satisfy the necessary number requirements, the ordinary differen- tial equation obtained in (2) will be a two point boundary value problem. One of the boundary conditions on the wall is known. Three of the remaining six are completely independent', and the other three are related in terms of constraint] conditions. The presence of . constraint] conditions makes the problem more difficult. In this re- port, the conventional initial value method was further developed, and a least square error method was proposed at the edge] of the boundary layer. It was found numerically that this method can be applied generally without regard for the number of unknown boundary conditions or the number of constraint| conditions. The problem of sublimating ablation was taken as an example to prove this. Calculations were performed by the above-described method, using Teflon as the ablation substance. The calculation conditions were the following: main istreamj Mach number 5-74, main stream] 38 stagnation point temperature 800° C - 1500° C, main ' Stream| stagna- tion point pressure 0.5 atm =2.0 atm, radius of curvature of the object 0.75 - 3.00 cm, and temperature inside the object of 20° C. The following conclusions were reached from the calculation results^ (i) Within the range of the : calculations performed here, the heat protection effect of ablation improves as the main jstreaml ' stagnation point temperature increases. This indicates the advan- tages of the ablation cooling method. (ii) When there is a greater radius of curvature of the object, the heat protection effect increases because the boundary layer be- comes thicker. (iii) The heat protection effect is better in solid models than in hollow models. However, thanks to the good thermal insula- tion properties of the ablation substance itself, the difference is not large, and there are possibilities that even substances with thin walls amounting to about 25^ of the radius of curvature can.-,-, be put into actual use. (iv) As for the ablation velocity distribution downstream in /22 the flow direction, there is a rather good coincidence when the main 'Stream] stagnation point temperature is above 1000° C, in comparison with the solutions obtained by the microdisturbance theory and with the test values. This indicates that this method can be applied effectively to problems of this type. It is expected that even better results could be attained if the two-dimensional effect of heat conduction inside the object were taken into consideration. However, in this case it would be neces- sary to use some method other than the H-W method, and the analysis at higher temperature ranges will ^be a task to be solved in the i future. 39 With reference to the solutions of the two point boundary value problems, the author wishes to express his heartfelt gratitude to director Horlkawa and specialist Watanabe of the Measuring Department for enabling the author to study the behavior of the solutions by means of an analog computer. APPENDIX A Boundary Layer Equations The following Is the boundary layer equation written In terms of the system of coordinates (X, Y) fixed In space for the flow pat- tern Illustrated In Figure 1. Mass conservation: d(pr.) ^ d {pUrc') 4. d (pV>.) ^Q dt SX dY Momentum conservation: 3< gx dY dX , dY V gy/ '' Ri, SY. (A.l) (A. 2) (A. 3) Energy conservation dt 3X dY 3j (A.il) 40 Diffusion: ^^-^^1|-^# (A. 5) Here, the quantities which have no suffixes 1, 2 are the mean quan- tities of a two-component mixed gas. Cp=C,iK+Cp2a-'K) (A. 6) The concentration is supposed to be K, e K. Since K-. + K^ = 1, it is true that Kp = 1 - K. Therefore, the analyses given below have been given a unified concentration of K. Since sweepback occurs on the surface of the object at a velo- city of V, J we can use the following transformation to coordinates (Xj y) fixed on the surface of the object: (A. 7) If we use the semi-stationary hypothesis, we obtain: dt Sy St 1^JL^--v,J- 9y (A. 8) By means of these, it is possible to transform Equations (A.l) through (A. 5) into a simultaneous equation at the system of coordi- nates (x, y ) . Furthermore, we can use the non-dlmensionallzed quantities expressed by: x=x/Rb, y=y/y*, r,=r,/Rt, u=u/u*, v=v/v*, p=p/p,i,, T=T/T,,H=H/Cp2T,, p='p/T.7u*, (A. 9) The- following non-dimensional simultaneous partial differential- equation can thus be obtained: 41 3 (pur,) y djpvr,) Bx By 3u Bu p„__+,„_^^ =0, Bx 3 / Bu \ ByK'' By ), Sp _ N (A. 10) (A. 11) (A. 12) 3H , „ BH pu— — +pv—-— Bx By B By BH Pr By 1 M 3"') ) -{(X-,),o<c„-c„>rf:} By BK , ,3K S / ^3K\ p=fiRT, H=C,T+ i Wtt\ (A. 13) (A.liJ) (A. 15) (A. 16) Let us assume that the shock wave layer is ;lnejllastlc and has a constant .density rj Since it is assumed that the velocity component on the surface of the object which is obtained from the inelastic hypersonic flow theory (Hayes and Probstein [21]_)[wlll be equal to u the velocity component at the end of the boundary layer. e' «,=«» V-y-yrfr''"'* (A. 17) Here (A. 18) Let us suppose that In this case. Furthermore, since '"^^^-- »,=sin*=a,(x)J H.=Cp.T.+ \Wu.^\ under conditions where H ' = .|1, r.=q^(i--tM'«.2)=r.w (A. 19) (A. 20) jA,_.2i)_ ■-■ (A. 22) ' 723 and the conditions at the boundary layer end will be given, 42 It is sufficient for the quantities to be non-dimensionalized as follows :| y* = (/<.» Rt/p.M*) ' V* =>* u*/St (A. 23) in order for the original boundary layer equation to be non-dimen- s ional and to assume the validity of Equations (A. 10) - (A.l4).| Using the state of the mean strean)],, let us rewrite the right side of y in the following manner ,r(2rAfco2-(r-l))^ kr-\) Moo^+i\KR,y (r+i) Mo,-,, i^=(;;«/p,»««/?.)'r-!^-— ■ -L_L:_^ Y V-^/(i+« • J (2rA/co2-(r-l)}^{(r-l)M«2+2}^ /?. "^ (r+l)MJ{8r/3(l+-i)2)^ (A. 24) V R. Furthermore , ifc, w=u*ycp2T,= y-^-/{i+x)Y ~<M«?aoc?ICp2T, (A. 25) (A. 26) In this case, we can use a,=ri?r,,c/.2/Cf2=r, c/»2-CD2=i?| to rewrite it as follows : (A. 27) 43 APPENDIX B Transfer Coefficient of Two-Component Gas Let us suppose that the first approximation of the Chapman- Enskog theory [22] applies to the two-component mixed gases which are being considered here. Let us then apply the Lennard-Jones (6, 12) potential to the Inter-partlcle forces. (1) Coefficient of viscosity /i=/<i/(l+Gi2-;i^ Ml 1-A" ^)+/<2/(l + G2,5J^^J3^)| (B.l) Here J (B.2) y-,-^an^ Up are the coefficients of viscosity of the component gases. expressed In the following terms: /i(=266. 93X 10-' (M,f ) * /ff,2j?,(2. 2>(gr.cm-i.sec-i) B.3) This Is a relationship which was derived for monoatomic gases, but it is known that it can be applied also to multlatomic gases as_w^J:.l_ [23]. Here, a. represents the radius of collision (A) and fi . ^^.'__n2__: the non-dimensional collision Integral. These are quantities which accompany the Lennard-Jones intermolecular function^ [ If we substi- tute Equation (B.3) in Equation (B.l) and convert it into non- dimensional terms, we obtain: + (l-z)/(l+(mG-l)z}](r/T.)5^ (B.4) Here, G = G21 , 11= m/fz, m = Mi/ Ml (B.5) 44 It Is assumed that the component gases present in the same space have the same temperature. (2) Thermal conduction coefficient 7: M. K Here , (B.6) (B.7) K-, and Kp are the thermal conduction coefficients for the component gases , /c, = -l^ -^ ft, (cal.cm-i.sec-'.oK-i) 4 Mi (B.:8) The Eucken factor Eu. Is a coefficient for approximating the thermal conduction coefficients of multlatomlc gases. It Is: Eu,=Q. 115+0. 2UCPiM,/R \ ( B . 9 ) When these are converted Into inon-dlmenslohal terms, we obtain: Here , + (l-z)/{l + (mS-l) ^}](Yr.)' (B.IO) /(0. 115* +0. 354M2C/>j) j 5=1.065G2, (B.ll) 45 (3) Chapman-Rubesln number /24 The Chapman-Rubesln number Is defined as : C=pil/p7it. \ (B.12) In view of Equation (3-9) and Equation (B.4), It can be expressed as follows : (b:13) 'I (4) Prandtl number Pr = ClflKr, Pr, Cp-!^/~ (B.14) Here, Prp Is the Prandtl number for the air at the external end of the boundary layer. In view of Equation (B.l4), Pr = Pr2 {].-o.z)\jz /{llG + a-fiG)z] + a-2)/(l + (mG-l)j] /[«/(pG+ (l-?5)zj+(l-2) /(n-(mG-])2)] (B.15) (5) Schmidt number Sc-'/ifpD =Sc.(:-(i-^)^}-£L(4.). (B.16) Here, SCp Is the Schmidt number of the air at the external end of the boundary layer. The two-component diffusion coefficient D Is (cm^.sec-') Ol2 = C<'l + 02)/2 (A), (B.17) 46 The non-dlmenslonal collision integrals were used for Q^ ^— In view of Equation (B.17), (B.18) + a-2)/{n-(mG-l)rj] References 1. Sutton, G. W. The Hydrodynamics and Heat Conduction of a Melting Surface. J. Aero/Space Sci., Vol. 25, 1958, pp. 29-32, 36. 2.- V Robert sVL-I Stagnation-Point Shielding by Melting and , ' Vaporlzation^l NASA TR R-10, 19 59- 3. Steg, L. and H. Lew. Hypersonic Ablation, The High Temp'erature Aspect of Hypersonic Plow. Pergamon Press, 1964, pp. 629-680 4. Swann, R. T. and J. South. A Theoretical Analysis of Effects of Ablation on Heat Transfer to an Arbitrary Axisymmetric Body. NASA TND-741, I96I. 5. Karashima, K., H. Kubota and K. Sato. An Aerodynamic Study of Ablation Near the Region of Stagnation Point of Axially Symmetric Bodies at Hypersonic Speeds,;] ISAS Rept . No. 425, Inst. Space and Aero. Sci., Univ.]Tokyo, I968. 6. Mills, A. F., A. V. Gomez and G. Strouhal. The Effect of Gas Phase Chemical Reactions on Heat Transfer to a Charring Ablator-;, AIAA Paper No. 70-869, 1970. 7. Nomura, Shlgeakl._. Abureshon ni yoru yodomi-ten kuriki kanetsu no gensho no'.sgkuitei (Measurement of Reduction of Stagnation Point Aerodynamic Heating by Means of Ablation). Koku Uchu Gijutsu Kenkyusho Hokoku (Technical Report of National Aerospace Laboratory), TR-I67, I968. Hartree, D. R. and J. R. Womersley. A Method for the Num'ericalf or Mechanical Solution of Certain Types of Partial Differen- tial Equations. Proc. Royal Soc . , I6IA, 1937, PP- 353-366. ~\ ;.l ^ 47 9. Smith, A. M. 0. and D. W. Clutter. Machine Calculation of Compressible Laminar Boundary Layers.] AIAA J., Vol. 3, No. 4, 1965, pp. 639-647. 10. Nachtheim, P. R. and P. Swlgert. Satisfaction of Asymptotic Boundary Conditions In Numerical Solution of Systems of Nonlinear Equations of Boundary-Layer Type. NASA TND-3004, 1965. 11. Fox J L., ed. Numerical Solution of Ordinary and Partial Differential Equations. Pergamon Press, 1962. 12. Lewallen, J. M. A Modified Quaslllnearlzatlon Concept for Solving the Nonlinear Two-Point Boundary Value Problem^ NASA TND-4025, 1967. 13. Bartlett, E. P., W. E. Nlcolet and J. T. Howe. Heat-Shield Ablation at Superorbltal Reentry Velocities j AIAA Paper 70-202, 1970. 14. Rashls, B. and R. N. Hopko. An Analytical Investigation of Ablation. NASA TMX-300 , I96O. /25 15. Kubota, H. An Aerodynamic Study of Hypersonic Heat Shield Problem with Local Mass Injection at Multiple Stations. ISAS Rept . No. 460, Inst. Space and Aero. Sci . , Univ. Tokyo, 1971. 16. Madorsky, S. L. Thermal Degradation of Organic Polymers. Polymer Reviews, Vol. 7, Intersclence, 1964, pp. 130-140. 17. Martin, J. J. Atmospheric Reentry^ An Introduction to Its Science and Engineering. Prentice-Hall, 1966. 18. Karashlma^ Keilchl, Kiyoshl Sato and Hlrotoshl Kubota. Kyokucho onsoku jiku-talsho donto buttal no yodoml-ten fukln nl okeru abureshon no jlkken (Ablation Experiments in the Vicinity of the Stagnation Point in Axlally Symmetrical Blunt-Headed Bodies at Hypersonic Speeds). Tokyo Daigaku Uchu Koku Kenkyusho Hokoku (Tokyo University, Institute of Space and Aeronautic Science, Reports), Vol. 4, No. 3 (A), 1968, pp. 335-347. 19. Halns , F. D. Equilibrium Shape of an Ablating Nose in Laminar Hypersonic Flow. AIAA J., Vol. 8, No. 7, 1970, pp. 1354- 1355. 20. Popper, L. A., T. Y. Toong and G. W. Sutton. Three-Dlmensional Ablation Considering Shape Changes and Internal Heat Conduc-- '' tlon. AIAA Paper No. 70-199, 1970; also Axlsymmetric Ablation with Shape Changes and Internal Heat -> ConductionTi AIAA J., Vol. 8, No. 11, 1970, pp. 2071-2074. 48 21. Hayes, W. D. and R. F. Probsteln. Hypersonic Flow Theory. Academic Press, 196O. ''22. Dorrance, W. H. Viscous Hypersonic Flow. McGraw Hill, 1962, 23. Bird, R. B., W. E. Stewart and E. N. Llghtfoot. Transport Phenomena. John Wiley & Sons, 19 60. Translated for National Aeronautics and Space Administration under c'ontract No. NASw 2035, by SCITRAN, P. 0. Box 5^56, Santa Barbara, California, 93108 49